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First-order transitions of system where both lattice site occupancy and lattice spacing fluctuate, 
such as cluster crystals, cannot be efficiently studied by traditional simulation methods, which 
necessarily fix one of these two degrees of freedom. The difficulty, however, can be surmounted 
by the generalized [N]pT ensemble [J. Chem. Phys. 136, 214106 (2012)]. Here we show that 
histogram reweighting and the [N]pT ensemble can be used to study an isostructural transition 
between cluster crystals of different occupancy in the generalized exponential model of index 4 
(GEM-4). Extending this scheme to finite-size scaling studies also allows to accurately determine 
the critical point parameters and to verify that it belongs to the Ising universality class. 
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Models with steeply growing but soft-core repulsions 
were first introduced as schematics for soft matter sys- 
tems, such as dendrimers and micellar crystals I n 
contrast to models with purely repulsive hard-core inter- 
actions, which simply crystallize, at high densities these 
models form cluster crystals in which each lattice site is 
multiply occupied. Clustering dramatically affects the 
materials properties of the crystal phase, contributing to 
both these model's response functions and elastic prop- 
erties . One of the canonical cluster-crystal formers, 
the generalized exponential model of index 4 (GEM-4) [f| 
with a pair interaction potential 

4>(r) =eexp[-(7» 4 ], (1) 

where e and a set the units of energy and length respec- 
tively, was also found to exhibit low-temperature first- 
order isostructural transitions between cluster phases 
with different integer occupancy, e.g., face-centered cubic 
(FCC) solids with double (FCC2) and triple (FCC3) oc- 
cupancy • Each of these isostructural transitions was 
found to terminate at a critical temperature T c , above 
which the lattice difference between the two phases van- 
ishes and particles are randomly distributed over the lat- 
tice sites. 

Although rather rare, isostructural transitions between 
singly occupied FCC solids with different lattice spac- 
ing were first observed in the cerium phase diagram over 
thirty years ago 0, HJ, and then in models with short- 
range attractive interactions HI QjJ. No direct study 
of their critical properties has been reported, but these 
properties are expected to be fairly standard, except in 
two dimensions, where coupling with a hexatic transi- 
tion may affect the critical behavior [llj. In the case 
of cluster crystals, however, the interplay between lat- 
tice occupancy and lattice spacing results in additional 
sources of fluctuation that could impact the critical be- 
havior, just like they affect these system's mechanical 
response. From a numerical simulation point of view, 
clustering also prevents the successful application of tra- 
ditional simulation methods. In this Brief Report, we 



present a finite-size scaling study of the FCC2-FCC3 crit- 
ical point in the GEM-4 by specially considering the lat- 
tice occupancy fluctuation at equilibrium. The approach 
confirms the Ising universality class of the transition and 
finds a minimal role for coupling. We expect that the ap- 
proach could be extended to other systems where lattice 
occupancy plays a central role, such as microphase form- 
ers [HI , crystals with large number of vacancies [H, [TiJ , 
and binary mixtures [TBI Il6|. 

Simulation techniques for locating first-order transi- 
tions and critical points (CP) are fairly well estab- 
lished [l8l - l2l1 ]. but there exists a key problem with sim- 
ulating cluster crystals. Once N c lattice sites are ini- 
tialized for a TV-particle system, the average lattice site 
occupancy n c = N/N c cannot generally relax to its equi- 
librium value. The recently developed [N]pT ensemble 
method successfully surmounts this problem by allowing 
both particle number and lattice spacing to fluctuate [IJ . 
At temperature T and pressure p, the equilibrium occu- 
pancy then coincides with the minimization of the con- 
strained Gibbs free energy G c (N,p, T, N c ) = fiN + fi c N c , 
whose exact differential form [22j ] 

dG c = -SdT + Vdp + ndN + f i c dN c , (2) 

depends not only on the entropy S, the volume V, the 
chemical potential of the particles /i, but also on /i c , the 
chemical potential-like quantity conjugate to N c . If N c 
is fixed, as it is in simulations of a reasonable size, the 
per particle free energy g c (p, T, n c ) = G c /N must thus be 
minimized to identify the equilibrium state at constant 
T and p. 

Coexistence densities. Given an external field Gn and 
a particle number window [iV m ; n , iV max ], the [N]pT en- 
semble samples N with probability 

V{N) ~ e^e~^, (3) 

and iterativcly determines g c within that window |4|. 
Convergence of the algorithm, which is identified by 
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FIG. 1: (Color online) Converged g c (black (thick) solid line) 
from [N]pT GEM-4 simulations with N c = 1372 at T = 0.045 
andp = 1.97. After reweighting to p — 1.9702 (red (thin) solid 
line), the two wells are at the same depth, which indicates that 
crystals at these two densities have the same chemical poten- 
tial /x eq (dotted line). Note that the g c has an intrinsic error 
of 0.0003, due to thermodynamic integration, which gives the 
pressure uncertainty of 0(0.001) and the optimized particle 
number uncertainty of 0(1). 



V(N) being flat, provides Gn differing from G c by a con- 
stant. This constant can then be determined by a single 
thermodynamic integration, neglecting the trivial ther- 
mal wavelength contribution [J, |23[ ■ Because the mini- 
mum of g c corresponds to the equilibrium lattice occu- 
pancy, with /i c — 0, the [N]pT ensemble method can 
also locate the coexistence densities of an isostructural 
transition by adjusting the temperature and pressure so 
that the two minima of g c have the same depth, i.e., the 
same p eq (Fig. [T|). If the simulated pressure differs from 
the coexistence pressure at the chosen temperature, his- 
togram reweighting over pressure to equate the two wells 
is used [H . By histogram reweighting over T and p, one 
can also obtain coexistence densities at neighboring phase 
points dHJ]. 

We first perform [N]pT simulations of a system with 
iV c = 1372 at T = 0.045 and p = 1.97. The pair po- 
tential 4>{r) is truncated at a cutoff distance of 1.7(7, be- 
yond which the potential energy is treated in an average 
way [23| . In addition to standard particle and logarithmic 
volume Monte Carlo (MC) moves [23], particle insertions 
(+) and removals (— ) are used with the acceptance ratio 
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± e /3AG ± -/3AE 



(4) 



where AG ± = G N ±i - G N , v + = V/(N + 1), r)~ = 
N/V, and AE^ is the energy cost of inserting/removing 
a particle [J, [25| . The free energy from thermodynamic 
integration is used as initial guess for Gn- For a particle 
window of width AN = 800, 10 7 MC cycles are sufficient 
to obtain well-averaged quantities. The resulting and 
reweighted g c curves are shown in Fig. Q] Although an 



intrinsic error in g c of the order 0(1/N) is introduced 
by the thermodynamic integration, the minimum of the 
g c curve is only affected by 5N ~ 0(1) [4|, which also 
sets the error bar on the coexistence pressure. With this 
scheme, we determine a series of coexistence densities 
for the FCC2-FCC3 isostructural transition close to the 
critical point (Fig. 2]). 

Critical point. The accurate determination of the crit- 
ical point involves a finite-size scaling approach based on 
the renormalization group |26j . At the apparent critical 
temperature T C (L) of a finite system of linear size L, the 
distribution of the order parameter adopts a universal 
form [2fjl |27|. once properly rescaled, that allows to ex- 
trapolate T c to the thermodynamic, i.e., infinite system 
size, limit. Note that under this definition, the distri- 
bution of order parameter still exhibits a double peak 
at the apparent critical temperature. In order to locate 
the GEM-4 FCC2-FCC3 isostructural critical point, we 
thus need to calculate the distribution V(p) of density 
p = N/V. We show below that one can obtain V(p) 
from G c — g c N by carefully considering the statistical 
mechanics of cluster crystals. 

At fixed N, p, and T (0 = l/k B T, where k B is the 
Boltzmann constant), the constrained Gibbs free energy 
G C {N C ) = G c {N 7 p,T 7 N c ) as a function of N c is effec- 
tively a Landau free energy with order parameter N c 



e -/3G c (JV c ) = ^ e -P G - 



- N c ), 



(5) 



where is a sum over all microstates indexed by v 
subject to the constant NpT constraint. The Kronecker 
delta function S(x) is unity when x = and zero oth- 
erwise. The equilibrium Gibbs free energy G(N,p,T) 
therefore satisfies 
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where the last approximation holds when the term with 
equilibrium A^ q dominates the sum. Analogously, we 
can define Landau free energy Gv(V) and G(N C , V) 



-l3Gv(V) _ 



^e-^S(V v -V), 



(7) 



and 



-0G(N c ,V) 



-PG, 



6(N» - N C )5{V V - V), (8) 



which arc related by 



-PG V (V) 



£« 



-PG(N C ,V) ^ pG{N^,V) 



(9) 



Here again, for a given volume V , there exists a specific 
NY , such that e~^ G ^ N " dominates the sum over N c . 
For a A^-particle system, density fluctuations result from 
changes in both V and N c . Yet because a single specific 
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FIG. 2: (Color online) The sampling distribution V(N,p) at 
T = 0.045 and p = 1.97 peaks at iVo for each density. On 
the p-N plane, the most sampled No (blue circles) is approx- 
imately equal to N = p(V) (red (thick) solid line). 



iVj dominates each V, we can consider that the density 
fluctuations are essentially along the N c = iVj contour 
of the two-dimensional G(N C , V) surface. In other words, 
the density fluctuation is governed by the Landau free 
energy G V (V) ~ G(N^ , V) for fixed NpT. 

In [NjpT ensemble simulations, N c is fixed while N 
fluctuates. The density distribution is then governed by 
a Landau free energy G(No(p), p) similarly defined as in 
Eq. [51 where Nq is the most probable particle number 
for a given density p. If we consider the [N]pT ensemble 
simulation to be a series of NpT simulations, then for 
each N E [A min , iV max ] there exists an average volume 
(V) . We can show that if densities are binned under the 
definition p = N/(V) for each N, then N « N for that 
density, with a relative error of the order 0{1/N). In the 
[N]pT ensemble, the conditional probability of observ- 
ing a particle number TV given the density p, V(N\p) ~ 
V(N,p) ~ NV(N, V), where the joint distribution 

V{N,V) ~ e PGAN.,p.T,N c ) e -p P V e -pF c {N.y,T,N c ) [| We 

can thus formally write 

e -/3(pV+F c (N,V,T,N )) e -l3G(N c ,V) 
V(N, V) e -/3Gc(N,p,T,N c ) = e -pd( Nc ,(V)) ' ^ 

Given pi, there exists N\ such that p\ — Ni/(V\). The 
relative probability of observing another N 2 = p 2 (V 2 } at 
pi is then 

V(N 2 , Pl ) _ N 2 V(N 2 ,N 2 / Pl ) _ e-^ N - N -IP^ 
V(N u pi) JVi P(N u {Vi)) ~ e -pG(N^v 2 )) < 

if the linear increase in probability due to the factor 
N 2 /N\ is negligible compared to the exponential suppres- 
sion. This condition is here obeyed because the standard 
deviation of N 2 from N\ is 0(1) (Fig. [2]). The most prob- 
able particle number N for a given density p — N/V is 
thus approximately the particle number under the def- 
inition N — p(V), within a relative error of the order 
0(1/N), as checked numerically in Fig. [51 We can there- 
fore use the Landau free energy G(N,p = N/(V)) — 



G C (N = p(V) ,p,T, N c ) to analyze the density fluctua- 
tions. 

Another challenge with [N]pT simulations is that both 
N and V are allowed to fluctuate, resulting in an ex- 
ponential growth of the probability distribution V(p = 

N/{V)) ~ e -PG c (N= P {v),p,T,N c )^ because f the i ncr ease 

in number of states with N. Although g c — G c /N shows 
a double well, the shape of G c is dominated by the linear 
increase with N, and the wells are nearly invisible on that 
scale. We can correct for this trivial exponential growth 
by normalizing the distribution with e - ^ ^ N , where p eq 
is the minimum of g c and thus the coexistence chemical 
potential. For a given system size L, the distribution 
of density is therefore P L (p = N/{V)) ~ ^(m^-Gc^ 
where G c = g c N is the converged field in the [N]pT sim- 
ulation. This approach is formally equivalent to running 
[N]pT simulation with Gn = p cq N without iterative up- 
dating (see Eq. [3]). This alternative strategy is, however, 
numerically inefficient if ^ eq is not known in advance or 
if the free energy barrier between the two minima is high 
and the sampling efficiency is low. This simulation ap- 
proach can thus be seen as one where the order parameter 
fluctuates at constant coexistence temperature, pressure 
and chemical potential. 

The linear system size L truncates the correlation 
length in a finite simulation box. In d-dimensional con- 
stant volume simulations, L is unambiguously defined as 
V x / d . When V fluctuates, such as in constant NpT sim- 
ulations, the fixed extensive quantity N is used as a mea- 
sure of length scale L oc iV 1 ^ [II], but in [NjpT simula- 
tions, neither N nor V are fixed. As mentioned earlier, 
the [N]pT simulation can be thought of as a series NpT 
simulations at various N's. One may thus be tempted 
to propose N 1 ^ as the linear size for each p = N/(V). 
Yet the resulting collapse is not good, because in cluster 
crystals the system size does not straightforwardly scale 
with N, due to clustering at fixed N c . We instead make 
an ansatz that L oc Nc . 

Finite-size scaling simulations. We perform [N]pT 
simulations for systems with iV c = 500, 864, andl372 at 
temperatures close to their T c (L)'s. For different sys- 

1 /3 

tern sizes L oc N c ' , the distributions Vl{p — N/(V)) ~ 
6 0(m N—Gc) a ^ their apparent critical temperature T C (L) 
as a function of the scaling variable, x — AL^/ U {p — 
p c {L)), collapse onto a universal function V(x), if the 
correct critical parameters T C (L), p c {L) and ft/v are cho- 
sen (26|. We use the Ising universality exponent 0/u = 
0.518 and the corresponding distribution V{x) [2!|. At 
T C (L), the ratio of the peak value of V(x) to its value at 
x = is about 2.173(4) [29]. With histogram reweight- 
ing, we identify the system-size dependent critical con- 
ditions (Table [T|, which result in good scaling behaviors 
(Fig. [3J). The nonuniversal normalization factor A is ar- 
bitrarily chosen as unity for the GEM-4 distributions, 
but other normalization conventions are also used in the 
literature (26j. Obtaining consistent data collapse using 
the empirical formula in Ref. [29| takes A = 0.132 for 
the Ising results. The infinite-system critical tempera- 
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FIG. 3: (Color online) Finite-size scaling of the GEM-4 
isostructural transition. The density distribution with ap- 
propriate critical values (Table [T| for the system of linear size 
L <x Nc /a , with N c = 500, 864, andl372, collapses onto the 
Ising universality function (J3/v = 0.518) [291 ]. Insets: the 
apparent critical temperature (left) and density (right) ex- 
trapolate to the infinite system size limit (Eq. I12|l . 
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FIG. 4: (Color online) GEM-4 first-order FCC2-FCC3 
isostructural transition. Previous thermodynamic integration 
(Tl) results [(| are shown for comparison. [N]pT simulations 
with histogram reweighting and finite-size scaling estimate 
T c (oo) = 0.0471(2) and p c (oo) = 1.2382(12). 



TABLE I: (Apparent) critical quantities for various system 
sizes lkJVc . 



N c 


T C (L) 


Pc(L) 


Pc (L) ^(L) 


500 


0.04555(10) 


1.2387(3) 


1.971(1) 2.9264(3) 


864 


0.0461(1) 


1.2386(3) 


1.972(1) 2.9298(3) 


1372 


0.0464(1) 


1.2385(3) 


1.973(1) 2.9316(3) 


oo 


0.0471(2) 


1.2382(12) 





ture and density can then be extracted from the scaling 
relations 

T c (L)-T c (oo)~L-^/» 

where the Ising universality class exponents are 9 = 0.54, 
a = 0.11 and v = 0.629 [lj| (Table I}. This analysis un- 
ambiguously assigns the GEM-4 isostructural transitions 
to the Ising universality class, and the accuracy of the 
critical point location is improved by an order of magni- 
tude (Fig. 13). 

In the liquid- vapor transition, which lacks the particle- 
hole symmetry observed in the Ising model, density p 
is not an appropriate order parameter to compare with 
the symmetric Ising model magnetization, so a linearly 
transformed density operator M = ^~ s " with parame- 
ters s and r, where u is the energy density, is needed [26| . 
For the GEM-4 isostructural transition, we find the sym- 
metry of p to be quite good close to T c . The relatively 
small density gap between the FCC 2 and FCC 3 phases 
(10-20%), compared to the orders of magnitude difference 
between vapor and liquid supports this observation. 

In summary, we connect the constrained Gibbs free en- 
ergy G c (N,p,T, N c ) with a Landau free energy G(N,p) 



to extract information about the density fluctuations of 
cluster crystals in the [N]pT ensemble. At fixed p, T 
and N c , the system can achieve a given density p by 
various combinations of N and V, corresponding to lat- 
tice occupancy and spacing fluctuations. Our analysis 
is based on the finding that a particular particle num- 
ber iV dominates at a given density and N can be ap- 
proximated by N = ply). If the distribution V(N,p) 
shown in Fig. [5] were to have broader peak over N, the 
lattice occupancy and thus density fluctuation would be- 
come stronger. In that case, one would need to perform 
random sampling in the [N]pT ensemble under a two- 
dimensional field G(N, V) to identify the minimum free 
energy contour, but the rest of the analysis would be sim- 
ilar as above. It would also possible to sample the lat- 
tice occupancy fluctuation by adapting the phase switch 
method [3l[ to the simulation of systems with fluctuating 
7V C [32J]. This approach, however, constrains particle oc- 
cupancy fluctuations - only planes of sites can be added 
or removed- which limits its accuracy, especially close to 
the critical point. 

The method and formalism above could also be applied 
to the study of binary mixtures of similar components, 
such as hard spheres whose diameter ratio cti/<T2 J$ 1.2. 
In their crystal form, these hard sphere mixtures exhibit 
first-order demixing between two FCC solids above a crit- 
ical pressure |15l - tL7| . Although liquid- liquid criticality 
can be studied using standard MC techniques [3(|, the 
fixed number of lattice sites in binary solids results in 
a problem similar to that of cluster crystals, where both 
the mole fraction of each component and the lattice spac- 
ing must allowed to fluctuate. Particle insertion/removal 
moves can be straightforwardly replaced by particle iden- 
tity changes under the Gibbs free energy per particle 
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g(x) = Gn/N = x/j,i+(l—x)[i2 as a function of mole frac- 
tion x, with the chemical potential of the two species [i\ 
and H2 mapped onto \i and \i c . Given T and p, the coex- 
istence is identified by the common tangent construction 
of the double-well g(x) curve. At the moment, the key 
difficulty in studying these systems is the low acceptance 
of growing small particles into large ones if the simula- 
tions are not sufficiently close to the critical point, but 
an improved initial guess for Gn could help solve this 
technical difficulty. 
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